Forecasting COVID-19 Cases and Deaths Using Time Series Analysis in R

# Run in console prior to running chunks:
getwd()
setwd('/Users/reginaduval/Grad_Work/MSDS692_Practicum1/Project/Data')
dist <- read.csv("County_data_distance.csv", header = TRUE, stringsAsFactors = FALSE)
train <- read.csv("train2.csv", header = TRUE, stringsAsFactors = FALSE)
pop <- read.csv("covid_county_population_usafacts.csv")
total_cases <- read.csv("total_cases.csv", header = TRUE)
total_deaths <- read.csv("total_deaths.csv", header = TRUE)
dist2 <- read.csv("Select_county_data_distance.csv", header = TRUE)
MN_9cum <- read.csv("MN_9_0613.csv", header = TRUE, stringsAsFactors = FALSE)
MN_9dcum <- read.csv("MN_9d_0613.csv", header = TRUE, stringsAsFactors = FALSE)
MN_0613 <- read.csv("covid_cases_0613.csv", header = TRUE, stringsAsFactors = FALSE)
county_data <- read.csv("County_data_distance2.csv", header = TRUE, stringsAsFactors = FALSE)
MN_0613d <- read.csv("covid_deaths_0613.csv", header = TRUE, stringsAsFactors = FALSE)
US_7cum <- read.csv("covid_cases_0613 copy.csv", header = TRUE, stringsAsFactors = FALSE)
US_7dcum <- read.csv("covid_deaths_0613 copy.csv", header = TRUE, stringsAsFactors = FALSE)
cases0612 <- read.csv("covid_confirmed_usafacts.csv", header = TRUE, stringsAsFactors = FALSE)
deaths0612 <- read.csv("covid_deaths_usafacts.csv", header = TRUE, stringsAsFactors = FALSE)

I. Data collection and clean up A. Pull initial data from Kaggle competition B. Investigate additional data points like population, distance, and meat-packing plants at county level C. Combine data and put into correct format for EDA and time series analysis

library(astsa)
library(broom)
library(caret)
library(DAAG)
library(dplyr)
library(dynlm)
library(forecast)
library(fpp2)
library(funModeling)
library(knitr)
library(MARSS)
library(matrixStats)
library(MTS)
library(quantmod)
library(readxl) 
library(tidyverse)
library(sf)
library(tidycensus)
library(tmap)
library(tmaptools)
library(tseries)
library(urbnmapr)
library(urca)
library(varhandle)
library(vars)

# # Review files from Kaggle competition
# train <- read.csv("train.csv", header = TRUE, stringsAsFactors = FALSE)
# # break out cases and deaths
# train_cases <- train %>% filter(Target == "ConfirmedCases")
# train_cases$Cases <- train_cases$TargetValue
# train_deaths <- train %>% filter(Target == "Fatalities")
# train_deaths$Deaths <- train_deaths$TargetValue
# # merge back to one dataset
# train2 <- merge(train_cases[, c(2, 3, 4, 5, 7, 10)], train_deaths[, c(2, 3, 4, 5, 7, 10)])
# # convert date column to date format (for EDA)
# train2$Date <- as.Date(train2$Date)
# # save train2
# write.csv(train2, "train2.csv", row.names = FALSE)

# # Create new file with extra variables (pop, dist, meat) by FIPS
# class(data)
# head(data)
# head(data2)
# tail(data3)
# data3$concatenate <- paste0(data3$county1, data3$county2)
# tail(data3)
# head(data4)
# data4 <- data4[, 1:6]
# head(data4)
# data4$concatenate <- data4$concatenate_BD
# head(data4)
# data4 <- data4[, -6]
# head(data4)
# dim(data4)
# data_full <- merge(data3[, c("concatenate", "mi_to_county")], data4, by.y, all.y = all)
# tail(data_full)
# summary(data_full)
# str(data_full)
# data_full$state <- as.character(data_full$state)
# data_full <- data_full %>% arrange(state)
# head(data_full)
# final_data <- data_full[, c(3, 4, 5, 6, 7, 2)]
# head(final_data)
# is.na(final_data)
# # write.csv(final_data, "County_data_distance.csv", row.names = FALSE)

# add pop to "County_data_distance.csv" from train2 dataset
# combine county data with train2 dataset
head(dist)

dim(dist)

[1] 3143   10

train_US <- train %>% filter(Country_Region == "US") %>% filter(County != "")
train_US$concatenate <- paste(train_US$County, train_US$Province_State, sep = "")
head(train_US)

dim(train_US)

[1] 339444      8

train_US <- train_US[, c(4, 8)]
dist2 <- left_join(dist, train_US, by = "concatenate")
dim(dist2)

[1] 337411     11

head(dist2)

dist2 <- unique(dist2)
dim(dist2)

[1] 3143   11

head(dist2)

# write.csv(dist2, "Distance Data/County_data_distance2.csv", row.names = FALSE)

2. Exploratory Data Analysis A. Basic EDA on county data B. Select specific counties to analyze more closely C. Add geographic data to map full US data and MN data

# COVID maps of US counties
tm_shape(spatial_data, projections = 2163) +
  tm_polygons(col = "total_cases", style = "quantile", palette = "OrRd") +
  tm_layout(title = "Total COVID cases by County",
            title.size = 1.1,
            title.position = c("center", "top"))
# COVID maps of US counties
tm_shape(spatial_data, projections = 2163) +
  tm_polygons(col = "total_cases", style = "quantile", palette = "OrRd") +
  tm_layout(title = "Total COVID cases by County",
            title.size = 1.1,
            title.position = c("center", "top"))

tm_shape(spatial_data, projections = 2163) +
  tm_polygons(col = "total_deaths", style = "log10_pretty", palette = "OrRd") +
  tm_layout(title = "Total COVID deaths by County",
            title.size = 1.1,
            title.position = c("center", "top"))

# Distance maps of MN counties
MNmap <- tm_shape(spatial_data_MN, projection = 2163) + 
  tm_polygons(col = "mi_to_county", style = "quantile", palette = "-YlGnBu") + 
  tm_legend(position = c("right", "center")) +
  tm_layout(title = "Minnesota Counties,\nDist to Twin Cities",
            title.size = 1.1,
            title.position = c("center", "top"))
# Distance maps of MN counties
MNmap <- tm_shape(spatial_data_MN, projection = 2163) + 
  tm_polygons(col = "mi_to_county", style = "quantile", palette = "-YlGnBu") + 
  tm_legend(position = c("right", "center")) +
  tm_layout(title = "Minnesota Counties,\nDist to Twin Cities",
            title.size = 1.1,
            title.position = c("center", "top"))
MNmap

MNmap + tm_shape(plants, projection = 2163) +
  tm_polygons(col = "meat_plant", palette = "Oranges")

# COVID maps of MN counties
MNmapCOVID <- tm_shape(spatial_data_MN, projection = 2163) + 
  tm_polygons(c("total_cases", "total_deaths"), style = c("quantile", "pretty"), 
              palette = list("YlGnBu", "YlOrRd")) + 
  tm_legend(position = c("right", "center")) +
  tm_layout(title = " MN Counties, COVID Data",
            title.size = 1.1,
            title.position = c("center", "top"))
MNmapCOVID

3. Country level time series analysis A. Filter data to 6 key countries (US, China, UK, Turkey, Brazil, India) B. Work through various time series models for US C. Choose best ts model (ARIMA) and tune D. Replicate models for other countries E. Forecast 10 days out

# BEST MODELS - use auto.arima models for simplicity and consistency
fit_casesBrazil <- auto.arima(casesBrazil, stepwise = FALSE, approximation = FALSE)
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_casesBrazil <- auto.arima(casesBrazil, stepwise = FALSE, approximation = FALSE)
fit_casesBrazil # ARIMA(3, 1, 2), AICc - 1680.4

Series: casesBrazil 
ARIMA(3,1,2) with drift 

Coefficients:
         ar1      ar2      ar3      ma1     ma2    drift
      0.2055  -0.1651  -0.6007  -0.5297  0.7085  96.2554
s.e.  0.1299   0.1944   0.0877   0.1086  0.2421  42.2660

sigma^2 estimated as 350911:  log likelihood=-832.63
AIC=1679.27   AICc=1680.4   BIC=1697.98

autoplot(fit_casesBrazil)

autoplot(fit_casesBrazil)

sarima.for(casesBrazil_train, n.ahead = 20, 3, 1, 2)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 6761.806 6826.082 6559.373 6232.534 6106.810 6311.566 6738.719 7131.645 7264.507 7098.194 6801.384
[12] 6627.818 6738.400 7092.859 7483.668 7684.478 7607.057 7358.058 7159.432 7191.605

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  383.3465  474.2080  562.6096  615.5042  654.6923  690.1179  732.5245  785.1342  842.0206  892.4678
[11]  930.7273  959.1391  984.2810 1012.8340 1048.6828 1089.9041 1129.6206 1162.0464 1186.8793 1208.0817

lines(casesBrazil_test)

checkresiduals(fit_casesBrazil)

    Ljung-Box test

data:  Residuals from ARIMA(3,1,2) with drift
Q* = 26.985, df = 16, p-value = 0.04165

Model df: 6.   Total lags used: 22

fit_casesBrazil2 <- arima(casesBrazil, order = c(4, 2, 4))
checkresiduals(fit_casesBrazil2) # passes the Ljung Box test

    Ljung-Box test

data:  Residuals from ARIMA(4,2,4)
Q* = 23.256, df = 14, p-value = 0.05622

Model df: 8.   Total lags used: 22

fit_deathsBrazil <- auto.arima(deathsBrazil, stepwise = FALSE, approximation = FALSE)
checkresiduals(fit_casesBrazil2) # passes the Ljung Box test

    Ljung-Box test

data:  Residuals from ARIMA(4,2,4)
Q* = 23.256, df = 14, p-value = 0.05622

Model df: 8.   Total lags used: 22

fit_deathsBrazil <- auto.arima(deathsBrazil, stepwise = FALSE, approximation = FALSE)
fit_deathsBrazil # ARIMA(4, 1, 0), AICc - 1155.41

Series: deathsBrazil 
ARIMA(4,1,0) with drift 

Coefficients:
          ar1      ar2      ar3      ar4   drift
      -0.2215  -0.0897  -0.1323  -0.4454  6.1765
s.e.   0.0914   0.1006   0.0985   0.0990  2.6035

sigma^2 estimated as 2642:  log likelihood=-571.28
AIC=1154.57   AICc=1155.41   BIC=1170.6

autoplot(fit_deathsBrazil)

autoplot(fit_deathsBrazil)

sarima.for(deathsBrazil_train, n.ahead = 20, 4, 1, 0)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 419.2291 347.3918 346.8175 394.1833 422.7578 436.1088 415.1895 399.8110 403.4813 419.4030 436.2125
[12] 440.1693 436.3332 433.4382 437.1307 446.1233 454.1783 458.2167 459.2964 460.8374

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 36.11087 44.27688 48.59298 49.47584 49.86034 52.18818 55.79384 59.64367 61.98798 63.32910 64.67596
[12] 66.46161 68.75936 70.98974 72.80179 74.31150 75.78037 77.40465 79.14916 80.85408

lines(deathsBrazil_test)

checkresiduals(fit_deathsBrazil) # passes

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0) with drift
Q* = 26.501, df = 17, p-value = 0.06581

Model df: 5.   Total lags used: 22

checkresiduals(fit_deathsBrazil) # passes

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0) with drift
Q* = 26.501, df = 17, p-value = 0.06581

Model df: 5.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_Brazil <- sarima.for(casesBrazil, n.ahead = 10, 4, 2, 4)

checkresiduals(fit_deathsBrazil) # passes

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0) with drift
Q* = 26.501, df = 17, p-value = 0.06581

Model df: 5.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_Brazil <- sarima.for(casesBrazil, n.ahead = 10, 4, 2, 4)

fc_10_Brazil$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 10178.22 10091.11 11590.62 12235.86 13126.08 13321.69 13724.37 14110.56 14851.42 15607.26

fcd_10_Brazil <- sarima.for(deathsBrazil, n.ahead = 10, 4, 1, 0) 

fcd_10_Brazil$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 643.2884 662.4851 594.1631 702.4151 691.7823 696.5816 724.2516 682.5469 705.0708 709.6913

checkresiduals(fit_casesBrazil2) # passes the Ljung Box test

    Ljung-Box test

data:  Residuals from ARIMA(4,2,4)
Q* = 23.256, df = 14, p-value = 0.05622

Model df: 8.   Total lags used: 22

fit_deathsBrazil <- auto.arima(deathsBrazil, stepwise = FALSE, approximation = FALSE)
fit_deathsBrazil # ARIMA(4, 1, 0), AICc - 1155.41

Series: deathsBrazil 
ARIMA(4,1,0) with drift 

Coefficients:
          ar1      ar2      ar3      ar4   drift
      -0.2215  -0.0897  -0.1323  -0.4454  6.1765
s.e.   0.0914   0.1006   0.0985   0.0990  2.6035

sigma^2 estimated as 2642:  log likelihood=-571.28
AIC=1154.57   AICc=1155.41   BIC=1170.6

autoplot(fit_deathsBrazil)

sarima.for(deathsBrazil_train, n.ahead = 20, 4, 1, 0)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 419.2291 347.3918 346.8175 394.1833 422.7578 436.1088 415.1895 399.8110 403.4813 419.4030 436.2125
[12] 440.1693 436.3332 433.4382 437.1307 446.1233 454.1783 458.2167 459.2964 460.8374

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 36.11087 44.27688 48.59298 49.47584 49.86034 52.18818 55.79384 59.64367 61.98798 63.32910 64.67596
[12] 66.46161 68.75936 70.98974 72.80179 74.31150 75.78037 77.40465 79.14916 80.85408

lines(deathsBrazil_test)

checkresiduals(fit_deathsBrazil) # passes

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0) with drift
Q* = 26.501, df = 17, p-value = 0.06581

Model df: 5.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_Brazil <- sarima.for(casesBrazil, n.ahead = 10, 4, 2, 4)

fc_10_Brazil$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 10178.22 10091.11 11590.62 12235.86 13126.08 13321.69 13724.37 14110.56 14851.42 15607.26

fcd_10_Brazil <- sarima.for(deathsBrazil, n.ahead = 10, 4, 1, 0) 

fcd_10_Brazil$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 643.2884 662.4851 594.1631 702.4151 691.7823 696.5816 724.2516 682.5469 705.0708 709.6913

# Repeat for China
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS, DIFFERENCED DATA SETS
casesChina <- tsChina[, "Cases"]
casesChina_train <- casesChina %>% window(end = c(2020, 120))
casesChina_train <- casesChina %>% window(end = c(2020, 120))
casesChina_test <- casesChina %>% window(end = c(2020, 130))
diff_casesChina <- diff(casesChina)
deathsChina <- tsChina[, "Deaths"]
deathsChina_train <- deathsChina %>% window(end = c(2020, 120))
deathsChina_test <- deathsChina %>% window(end = c(2020, 130))
diff_deathsChina <- diff(deathsChina)
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_casesChina <- auto.arima(casesChina, stepwise = FALSE, approximation = FALSE)
fit_casesChina # ARIMA(0, 1, 1), AICc - 1874.66

Series: casesChina 
ARIMA(0,1,1) 

Coefficients:
          ma1
      -0.7238
s.e.   0.0734

sigma^2 estimated as 2295057:  log likelihood=-935.27
AIC=1874.54   AICc=1874.66   BIC=1879.89

fit_casesChina # ARIMA(0, 1, 1), AICc - 1874.66

Series: casesChina 
ARIMA(0,1,1) 

Coefficients:
          ma1
      -0.7238
s.e.   0.0734

sigma^2 estimated as 2295057:  log likelihood=-935.27
AIC=1874.54   AICc=1874.66   BIC=1879.89

autoplot(fit_casesChina)

autoplot(fit_casesChina)

sarima.for(casesChina_train, n.ahead = 20, 0, 1, 1)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  -14.44721  -22.21844  -29.98968  -37.76091  -45.53214  -53.30338  -61.07461  -68.84585  -76.61708
[10]  -84.38831  -92.15955  -99.93078 -107.70202 -115.47325 -123.24449 -131.01572 -138.78695 -146.55819
[19] -154.32942 -162.10066

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 1583.415 1642.436 1699.409 1754.533 1807.976 1859.885 1910.384 1959.582 2007.574 2054.446 2100.272
[12] 2145.119 2189.048 2232.112 2274.361 2315.839 2356.588 2396.644 2436.041 2474.811

lines(casesChina_test)

checkresiduals(fit_casesChina) # passes Ljung Box test

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)
Q* = 10.617, df = 21, p-value = 0.9697

Model df: 1.   Total lags used: 22

fit_deathsChina <- auto.arima(deathsChina, stepwise = FALSE, approximation = FALSE)
fit_deathsChina # ARIMA(0, 0, 0) WHITE NOISE,AICc - 1361.29

Series: deathsChina 
ARIMA(0,0,0) with non-zero mean 

Coefficients:
         mean
      42.7778
s.e.  12.4686

sigma^2 estimated as 16947:  log likelihood=-678.59
AIC=1361.18   AICc=1361.29   BIC=1366.54

# Repeat for India
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS, DIFFERENCED DATA SETS
casesIndia <- tsIndia[, "Cases"]
casesIndia_train <- casesIndia %>% window(end = c(2020, 120))
# Repeat for India
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS, DIFFERENCED DATA SETS
casesIndia <- tsIndia[, "Cases"]
casesIndia_train <- casesIndia %>% window(end = c(2020, 120))
casesIndia_test <- casesIndia %>% window(end = c(2020, 130))
diff_casesIndia <- diff(casesIndia)
diff_casesIndia <- diff(casesIndia)
deathsIndia <- tsIndia[, "Deaths"]
deathsIndia <- tsIndia[, "Deaths"]
deathsIndia_train <- deathsIndia %>% window(end = c(2020, 120))
deathsIndia_train <- deathsIndia %>% window(end = c(2020, 120))
deathsIndia_test <- deathsIndia %>% window(end = c(2020, 130))
diff_deathsIndia <- diff(deathsIndia)
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_casesIndia <- auto.arima(casesIndia, stepwise = FALSE, approximation = FALSE)
fit_casesIndia # ARIMA(0, 1, 2), AICc - 1470.79

Series: casesIndia 
ARIMA(0,1,2) with drift 

Coefficients:
          ma1     ma2    drift
      -0.6010  0.2519  29.7706
s.e.   0.0972  0.0920  14.1403

sigma^2 estimated as 51751:  log likelihood=-731.2
AIC=1470.4   AICc=1470.79   BIC=1481.09

autoplot(fit_casesIndia)

sarima.for(casesIndia_train, n.ahead = 20, 0, 1, 2)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 1688.654 1840.992 1859.485 1877.977 1896.470 1914.963 1933.456 1951.949 1970.441 1988.934 2007.427
[12] 2025.920 2044.413 2062.905 2081.398 2099.891 2118.384 2136.877 2155.369 2173.862

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 155.6078 173.5986 181.8841 189.8083 197.4147 204.7387 211.8096 218.6520 225.2866 231.7314 238.0017
[12] 244.1110 250.0711 255.8924 261.5842 267.1547 272.6115 277.9611 283.2098 288.3629

lines(casesIndia_test)

checkresiduals(fit_casesIndia) # passes the Ljung Box test

    Ljung-Box test

data:  Residuals from ARIMA(0,1,2) with drift
Q* = 23.363, df = 19, p-value = 0.2217

Model df: 3.   Total lags used: 22

fit_deathsIndia <- auto.arima(deathsIndia, stepwise = FALSE, approximation = FALSE)
fit_deathsIndia <- auto.arima(deathsIndia, stepwise = FALSE, approximation = FALSE)
fit_deathsIndia # ARIMA(4, 1, 0), AICc - 848.95

Series: deathsIndia 
ARIMA(4,1,0) with drift 

Coefficients:
          ar1      ar2      ar3      ar4   drift
      -0.7113  -0.4206  -0.2434  -0.2459  1.0281
s.e.   0.0940   0.1145   0.1139   0.0955  0.4478

sigma^2 estimated as 151:  log likelihood=-418.06
AIC=848.11   AICc=848.95   BIC=864.15

fit_deathsIndia # ARIMA(4, 1, 0), AICc - 848.95

Series: deathsIndia 
ARIMA(4,1,0) with drift 

Coefficients:
          ar1      ar2      ar3      ar4   drift
      -0.7113  -0.4206  -0.2434  -0.2459  1.0281
s.e.   0.0940   0.1145   0.1139   0.0955  0.4478

sigma^2 estimated as 151:  log likelihood=-418.06
AIC=848.11   AICc=848.95   BIC=864.15

autoplot(fit_deathsIndia)

autoplot(fit_deathsIndia)

sarima.for(deathsIndia_train, n.ahead = 20, 4, 1, 0)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 62.68084 67.86632 66.96706 67.84111 70.50516 68.72801 70.70753 71.23923 71.34247 72.92316 72.91385
[12] 73.83024 74.66868 75.02204 76.01576 76.54349 77.21679 78.01620 78.57044 79.33989

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  6.254629  6.371803  6.690971  7.295919  7.393912  8.038123  8.313837  8.559386  8.974898  9.180993
[11]  9.505059  9.788530 10.020246 10.315512 10.548093 10.797000 11.048356 11.271317 11.511939 11.735595

lines(deathsIndia_test)

lines(deathsIndia_test)

checkresiduals(fit_deathsIndia) # passes

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0) with drift
Q* = 9.9658, df = 17, p-value = 0.905

Model df: 5.   Total lags used: 22

checkresiduals(fit_deathsIndia) # passes

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0) with drift
Q* = 9.9658, df = 17, p-value = 0.905

Model df: 5.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_India <- sarima.for(casesIndia, n.ahead = 10, 0, 1, 2)

# Use best models to forecast further ahead
fc_10_India <- sarima.for(casesIndia, n.ahead = 10, 0, 1, 2)

fc_10_India$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 3296.262 3246.772 3276.542 3306.313 3336.083 3365.854 3395.624 3425.395 3455.166 3484.936

fcd_10_India <- sarima.for(deathsIndia, n.ahead = 10, 4, 1, 0) 

fcd_10_India$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 113.5211 108.5624 112.9270 110.2884 114.8406 115.5648 115.3987 117.4482 117.4596 119.1467

# Repeat for Turkey
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS, DIFFERENCED DATA SETS
casesTurkey <- tsTurkey[, "Cases"]
casesTurkey_train <- casesTurkey %>% window(end = c(2020, 120))
# Repeat for Turkey
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS, DIFFERENCED DATA SETS
casesTurkey <- tsTurkey[, "Cases"]
casesTurkey_train <- casesTurkey %>% window(end = c(2020, 120))
casesTurkey_test <- casesTurkey %>% window(end = c(2020, 130))
diff_casesTurkey <- diff(casesTurkey)
diff_casesTurkey <- diff(casesTurkey)
deathsTurkey <- tsTurkey[, "Deaths"]
deathsTurkey <- tsTurkey[, "Deaths"]
deathsTurkey_train <- deathsTurkey %>% window(end = c(2020, 120))
fcd_10_India$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 113.5211 108.5624 112.9270 110.2884 114.8406 115.5648 115.3987 117.4482 117.4596 119.1467

# Repeat for Turkey
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS, DIFFERENCED DATA SETS
casesTurkey <- tsTurkey[, "Cases"]
casesTurkey_train <- casesTurkey %>% window(end = c(2020, 120))
casesTurkey_test <- casesTurkey %>% window(end = c(2020, 130))
diff_casesTurkey <- diff(casesTurkey)
deathsTurkey <- tsTurkey[, "Deaths"]
deathsTurkey_train <- deathsTurkey %>% window(end = c(2020, 120))
deathsTurkey_test <- deathsTurkey %>% window(end = c(2020, 130))
casesTurkey_test <- casesTurkey %>% window(end = c(2020, 130))
diff_casesTurkey <- diff(casesTurkey)
deathsTurkey <- tsTurkey[, "Deaths"]
deathsTurkey_train <- deathsTurkey %>% window(end = c(2020, 120))
deathsTurkey_test <- deathsTurkey %>% window(end = c(2020, 130))
diff_deathsTurkey <- diff(deathsTurkey)
diff_deathsTurkey <- diff(deathsTurkey)
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_casesTurkey <- auto.arima(casesTurkey, stepwise = FALSE, approximation = FALSE)
fit_casesTurkey # ARIMA(1, 1, 4), AICc - 1518.22

Series: casesTurkey 
ARIMA(1,1,4) 

Coefficients:
         ar1      ma1      ma2      ma3     ma4
      0.8292  -0.8253  -0.1877  -0.1799  0.5037
s.e.  0.0804   0.1027   0.1077   0.1282  0.0977

sigma^2 estimated as 77758:  log likelihood=-752.69
AIC=1517.38   AICc=1518.22   BIC=1533.42

autoplot(fit_casesTurkey)

autoplot(fit_casesTurkey)

sarima.for(casesTurkey_train, n.ahead = 20, 1, 1, 4)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 3029.565 3008.048 2964.522 3387.826 3483.513 3527.171 3562.567 3596.650 3630.524 3664.366 3698.202
[12] 3732.038 3765.873 3799.709 3833.544 3867.380 3901.215 3935.050 3968.886 4002.721

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  279.6939  400.2613  503.6633  560.9737  677.8426  787.9286  885.9506  974.3812 1055.4610 1130.7467
[11] 1201.3243 1267.9796 1331.3018 1391.7459 1449.6720 1505.3708 1559.0810 1611.0015 1661.3001 1710.1199

lines(casesTurkey_test)

checkresiduals(fit_casesTurkey) # passes Ljung Box test

    Ljung-Box test

data:  Residuals from ARIMA(1,1,4)
Q* = 24.078, df = 17, p-value = 0.1173

Model df: 5.   Total lags used: 22

checkresiduals(fit_casesTurkey) # passes Ljung Box test

    Ljung-Box test

data:  Residuals from ARIMA(1,1,4)
Q* = 24.078, df = 17, p-value = 0.1173

Model df: 5.   Total lags used: 22

fit_deathsTurkey <- auto.arima(deathsTurkey, stepwise = FALSE, approximation = FALSE)
fit_deathsTurkey <- auto.arima(deathsTurkey, stepwise = FALSE, approximation = FALSE)
fit_deathsTurkey # ARIMA(1, 1, 1), AICc - 648.99

Series: deathsTurkey 
ARIMA(1,1,1) 

Coefficients:
         ar1      ma1
      0.9280  -0.7925
s.e.  0.0567   0.0843

sigma^2 estimated as 24.17:  log likelihood=-321.38
AIC=648.76   AICc=648.99   BIC=656.77

fit_deathsTurkey # ARIMA(1, 1, 1), AICc - 648.99

Series: deathsTurkey 
ARIMA(1,1,1) 

Coefficients:
         ar1      ma1
      0.9280  -0.7925
s.e.  0.0567   0.0843

sigma^2 estimated as 24.17:  log likelihood=-321.38
AIC=648.76   AICc=648.99   BIC=656.77

autoplot(fit_deathsTurkey)

sarima.for(deathsTurkey_train, n.ahead = 20, 1, 1, 3)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  88.75692  88.90431  89.39671  90.28286  91.17014  92.05743  92.94471  93.83200  94.71928  95.60657
[11]  96.49385  97.38114  98.26842  99.15571 100.04299 100.93028 101.81756 102.70484 103.59213 104.47941

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  4.500711  6.979740  9.259851 11.696584 13.708356 15.460534 17.033414 18.472852 19.807962 21.058595
[11] 22.239009 23.359850 24.429320 25.453894 26.438794 27.388299 28.305971 29.194812 30.057381 30.895877

lines(deathsTurkey_test)

checkresiduals(fit_deathsTurkey) # passes

    Ljung-Box test

data:  Residuals from ARIMA(1,1,1)
Q* = 23.582, df = 20, p-value = 0.2611

Model df: 2.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_Turkey <- sarima.for(casesTurkey, n.ahead = 10, 1, 1, 4)

fc_10_Turkey$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 1670.783 1717.924 1816.134 1763.657 1722.394 1690.408 1666.100 1648.142 1635.438 1627.083

fcd_10_Turkey <- sarima.for(deathsTurkey, n.ahead = 10, 1, 1, 3) 

fcd_10_Turkey$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 47.58606 45.80776 44.06012 42.45847 40.99171 39.64960 38.42267 37.30217 36.28003 35.34877

# Repeat for UK
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS, DIFFERENCED DATA SETS
casesUK <- tsUK[, "Cases"]
casesUK_train <- casesUK %>% window(end = c(2020, 120))
casesUK_test <- casesUK %>% window(end = c(2020, 130))
casesUK_test <- casesUK %>% window(end = c(2020, 130))
diff_casesUK <- diff(casesUK)
deathsUK <- tsUK[, "Deaths"]
sarima.for(deathsTurkey_train, n.ahead = 20, 1, 1, 3)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  88.75692  88.90431  89.39671  90.28286  91.17014  92.05743  92.94471  93.83200  94.71928  95.60657
[11]  96.49385  97.38114  98.26842  99.15571 100.04299 100.93028 101.81756 102.70484 103.59213 104.47941

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  4.500711  6.979740  9.259851 11.696584 13.708356 15.460534 17.033414 18.472852 19.807962 21.058595
[11] 22.239009 23.359850 24.429320 25.453894 26.438794 27.388299 28.305971 29.194812 30.057381 30.895877

lines(deathsTurkey_test)

checkresiduals(fit_deathsTurkey) # passes

    Ljung-Box test

data:  Residuals from ARIMA(1,1,1)
Q* = 23.582, df = 20, p-value = 0.2611

Model df: 2.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_Turkey <- sarima.for(casesTurkey, n.ahead = 10, 1, 1, 4)

fc_10_Turkey$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 1670.783 1717.924 1816.134 1763.657 1722.394 1690.408 1666.100 1648.142 1635.438 1627.083

fcd_10_Turkey <- sarima.for(deathsTurkey, n.ahead = 10, 1, 1, 3) 

fcd_10_Turkey$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 47.58606 45.80776 44.06012 42.45847 40.99171 39.64960 38.42267 37.30217 36.28003 35.34877

# Repeat for UK
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS, DIFFERENCED DATA SETS
casesUK <- tsUK[, "Cases"]
casesUK_train <- casesUK %>% window(end = c(2020, 120))
casesUK_test <- casesUK %>% window(end = c(2020, 130))
diff_casesUK <- diff(casesUK)
deathsUK <- tsUK[, "Deaths"]
deathsUK_train <- deathsUK %>% window(end = c(2020, 120))
diff_casesUK <- diff(casesUK)
deathsUK <- tsUK[, "Deaths"]
deathsUK_train <- deathsUK %>% window(end = c(2020, 120))
deathsUK_test <- deathsUK %>% window(end = c(2020, 130))
deathsUK_test <- deathsUK %>% window(end = c(2020, 130))
diff_deathsUK <- diff(deathsUK)
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_casesUK <- auto.arima(casesUK, stepwise = FALSE, approximation = FALSE)
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_casesUK <- auto.arima(casesUK, stepwise = FALSE, approximation = FALSE)
fit_casesUK # ARIMA(3, 1, 2), AICc - 1700.4

Series: casesUK 
ARIMA(3,1,2) with drift 

Coefficients:
         ar1      ar2      ar3      ma1     ma2    drift
      0.6527  -0.2495  -0.3920  -1.3253  0.7142  43.0059
s.e.  0.1770   0.1084   0.1273   0.1945  0.1163  24.2431

sigma^2 estimated as 423776:  log likelihood=-842.63
AIC=1699.27   AICc=1700.4   BIC=1717.98

fit_casesUK # ARIMA(3, 1, 2), AICc - 1700.4

Series: casesUK 
ARIMA(3,1,2) with drift 

Coefficients:
         ar1      ar2      ar3      ma1     ma2    drift
      0.6527  -0.2495  -0.3920  -1.3253  0.7142  43.0059
s.e.  0.1770   0.1084   0.1273   0.1945  0.1163  24.2431

sigma^2 estimated as 423776:  log likelihood=-842.63
AIC=1699.27   AICc=1700.4   BIC=1717.98

autoplot(fit_casesUK)

fit_casesUK # ARIMA(3, 1, 2), AICc - 1700.4

Series: casesUK 
ARIMA(3,1,2) with drift 

Coefficients:
         ar1      ar2      ar3      ma1     ma2    drift
      0.6527  -0.2495  -0.3920  -1.3253  0.7142  43.0059
s.e.  0.1770   0.1084   0.1273   0.1945  0.1163  24.2431

sigma^2 estimated as 423776:  log likelihood=-842.63
AIC=1699.27   AICc=1700.4   BIC=1717.98

autoplot(fit_casesUK)

sarima.for(casesUK_train, n.ahead = 20, 3, 1, 2)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 4050.904 4152.183 4231.449 4326.202 4377.655 4403.729 4408.403 4419.534 4449.988 4504.446 4572.139
[12] 4638.233 4690.415 4725.934 4751.002 4776.431 4811.093 4857.667 4912.160 4967.069

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  586.0736  606.6333  662.6996  676.5956  711.1349  753.3029  818.8537  884.9817  942.5691  983.4384
[11] 1012.5957 1036.5155 1061.9516 1093.0877 1130.5595 1170.6427 1208.1820 1240.1062 1266.7283 1290.4919

lines(casesUK_test)

lines(casesUK_test)

checkresiduals(fit_casesUK) # passes the LJung Box test

    Ljung-Box test

data:  Residuals from ARIMA(3,1,2) with drift
Q* = 22.681, df = 16, p-value = 0.1225

Model df: 6.   Total lags used: 22

checkresiduals(fit_casesUK) # passes the LJung Box test

    Ljung-Box test

data:  Residuals from ARIMA(3,1,2) with drift
Q* = 22.681, df = 16, p-value = 0.1225

Model df: 6.   Total lags used: 22

fit_deathsUK <- auto.arima(deathsUK, stepwise = FALSE, approximation = FALSE)
fit_deathsUK # ARIMA(1, 1, 3), AICc - 1350.23

Series: deathsUK 
ARIMA(1,1,3) 

Coefficients:
         ar1      ma1      ma2     ma3
      0.5772  -1.1140  -0.1400  0.5424
s.e.  0.1186   0.1062   0.1455  0.0839

sigma^2 estimated as 16286:  log likelihood=-669.82
AIC=1349.64   AICc=1350.23   BIC=1363

autoplot(fit_deathsUK)

autoplot(fit_deathsUK)

sarima.for(deathsUK_train, n.ahead = 20, 1, 1, 3)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 595.2481 585.6878 551.3371 532.9061 523.9270 520.5598 520.5246 522.4676 525.5851 529.4001 533.6291
[12] 538.1039 542.7246 547.4320 552.1909 556.9803 561.7879 566.6062 571.4309 576.2594

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 123.3315 136.4543 136.7248 142.5479 154.3391 169.6170 186.2735 203.0971 219.5073 235.2700 250.3215
[12] 264.6768 278.3825 291.4956 304.0729 316.1669 327.8245 339.0871 349.9906 360.5663

autoplot(fit_deathsUK)

sarima.for(deathsUK_train, n.ahead = 20, 1, 1, 3)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 595.2481 585.6878 551.3371 532.9061 523.9270 520.5598 520.5246 522.4676 525.5851 529.4001 533.6291
[12] 538.1039 542.7246 547.4320 552.1909 556.9803 561.7879 566.6062 571.4309 576.2594

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 123.3315 136.4543 136.7248 142.5479 154.3391 169.6170 186.2735 203.0971 219.5073 235.2700 250.3215
[12] 264.6768 278.3825 291.4956 304.0729 316.1669 327.8245 339.0871 349.9906 360.5663

lines(deathsUK_test)

checkresiduals(fit_deathsUK)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,3)
Q* = 70.983, df = 18, p-value = 3.079e-08

Model df: 4.   Total lags used: 22

checkresiduals(fit_deathsUK)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,3)
Q* = 70.983, df = 18, p-value = 3.079e-08

Model df: 4.   Total lags used: 22

fit_deathsUK2 <- arima(deathsUK, order = c(6, 2, 9))
checkresiduals(fit_deathsUK)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,3)
Q* = 70.983, df = 18, p-value = 3.079e-08

Model df: 4.   Total lags used: 22

fit_deathsUK2 <- arima(deathsUK, order = c(6, 2, 9))
checkresiduals(fit_deathsUK2) # p-value = 0.019, closest I could get to 0.05

    Ljung-Box test

data:  Residuals from ARIMA(6,2,9)
Q* = 16.686, df = 7, p-value = 0.01953

Model df: 15.   Total lags used: 22

fit_deathsUK2 <- arima(deathsUK, order = c(6, 2, 9))
checkresiduals(fit_deathsUK2) # p-value = 0.019, closest I could get to 0.05

    Ljung-Box test

data:  Residuals from ARIMA(6,2,9)
Q* = 16.686, df = 7, p-value = 0.01953

Model df: 15.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_UK <- sarima.for(casesUK, n.ahead = 10, 3, 1, 2)

fc_10_UK$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 4248.787 4537.199 4975.123 5093.243 4990.550 4764.910 4639.472 4696.669 4896.270 5103.977

fcd_10_UK <- sarima.for(deathsUK, n.ahead = 10, 6, 2, 9) 

fit_deathsTurkey <- auto.arima(deathsTurkey, stepwise = FALSE, approximation = FALSE)
fit_deathsTurkey # ARIMA(1, 1, 1), AICc - 648.99

Series: deathsTurkey 
ARIMA(1,1,1) 

Coefficients:
         ar1      ma1
      0.9280  -0.7925
s.e.  0.0567   0.0843

sigma^2 estimated as 24.17:  log likelihood=-321.38
AIC=648.76   AICc=648.99   BIC=656.77

autoplot(fit_deathsTurkey)

sarima.for(deathsTurkey_train, n.ahead = 20, 1, 1, 3)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  88.75692  88.90431  89.39671  90.28286  91.17014  92.05743  92.94471  93.83200  94.71928  95.60657
[11]  96.49385  97.38114  98.26842  99.15571 100.04299 100.93028 101.81756 102.70484 103.59213 104.47941

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  4.500711  6.979740  9.259851 11.696584 13.708356 15.460534 17.033414 18.472852 19.807962 21.058595
[11] 22.239009 23.359850 24.429320 25.453894 26.438794 27.388299 28.305971 29.194812 30.057381 30.895877

lines(deathsTurkey_test)

checkresiduals(fit_deathsTurkey) # passes

    Ljung-Box test

data:  Residuals from ARIMA(1,1,1)
Q* = 23.582, df = 20, p-value = 0.2611

Model df: 2.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_Turkey <- sarima.for(casesTurkey, n.ahead = 10, 1, 1, 4)

fc_10_Turkey$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 1670.783 1717.924 1816.134 1763.657 1722.394 1690.408 1666.100 1648.142 1635.438 1627.083

fcd_10_Turkey <- sarima.for(deathsTurkey, n.ahead = 10, 1, 1, 3) 

fcd_10_Turkey$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 47.58606 45.80776 44.06012 42.45847 40.99171 39.64960 38.42267 37.30217 36.28003 35.34877

# Repeat for UK
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS, DIFFERENCED DATA SETS
casesUK <- tsUK[, "Cases"]
casesUK_train <- casesUK %>% window(end = c(2020, 120))
casesUK_test <- casesUK %>% window(end = c(2020, 130))
diff_casesUK <- diff(casesUK)
deathsUK <- tsUK[, "Deaths"]
deathsUK_train <- deathsUK %>% window(end = c(2020, 120))
deathsUK_test <- deathsUK %>% window(end = c(2020, 130))
diff_deathsUK <- diff(deathsUK)
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_casesUK <- auto.arima(casesUK, stepwise = FALSE, approximation = FALSE)
fit_casesUK # ARIMA(3, 1, 2), AICc - 1700.4

Series: casesUK 
ARIMA(3,1,2) with drift 

Coefficients:
         ar1      ar2      ar3      ma1     ma2    drift
      0.6527  -0.2495  -0.3920  -1.3253  0.7142  43.0059
s.e.  0.1770   0.1084   0.1273   0.1945  0.1163  24.2431

sigma^2 estimated as 423776:  log likelihood=-842.63
AIC=1699.27   AICc=1700.4   BIC=1717.98

autoplot(fit_casesUK)

sarima.for(casesUK_train, n.ahead = 20, 3, 1, 2)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 4050.904 4152.183 4231.449 4326.202 4377.655 4403.729 4408.403 4419.534 4449.988 4504.446 4572.139
[12] 4638.233 4690.415 4725.934 4751.002 4776.431 4811.093 4857.667 4912.160 4967.069

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  586.0736  606.6333  662.6996  676.5956  711.1349  753.3029  818.8537  884.9817  942.5691  983.4384
[11] 1012.5957 1036.5155 1061.9516 1093.0877 1130.5595 1170.6427 1208.1820 1240.1062 1266.7283 1290.4919

lines(casesUK_test)

checkresiduals(fit_casesUK) # passes the LJung Box test

    Ljung-Box test

data:  Residuals from ARIMA(3,1,2) with drift
Q* = 22.681, df = 16, p-value = 0.1225

Model df: 6.   Total lags used: 22

fit_deathsUK <- auto.arima(deathsUK, stepwise = FALSE, approximation = FALSE)
fit_deathsUK # ARIMA(1, 1, 3), AICc - 1350.23

Series: deathsUK 
ARIMA(1,1,3) 

Coefficients:
         ar1      ma1      ma2     ma3
      0.5772  -1.1140  -0.1400  0.5424
s.e.  0.1186   0.1062   0.1455  0.0839

sigma^2 estimated as 16286:  log likelihood=-669.82
AIC=1349.64   AICc=1350.23   BIC=1363

autoplot(fit_deathsUK)

sarima.for(deathsUK_train, n.ahead = 20, 1, 1, 3)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 595.2481 585.6878 551.3371 532.9061 523.9270 520.5598 520.5246 522.4676 525.5851 529.4001 533.6291
[12] 538.1039 542.7246 547.4320 552.1909 556.9803 561.7879 566.6062 571.4309 576.2594

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 123.3315 136.4543 136.7248 142.5479 154.3391 169.6170 186.2735 203.0971 219.5073 235.2700 250.3215
[12] 264.6768 278.3825 291.4956 304.0729 316.1669 327.8245 339.0871 349.9906 360.5663

lines(deathsUK_test)

checkresiduals(fit_deathsUK)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,3)
Q* = 70.983, df = 18, p-value = 3.079e-08

Model df: 4.   Total lags used: 22

fit_deathsUK2 <- arima(deathsUK, order = c(6, 2, 9))
checkresiduals(fit_deathsUK2) # p-value = 0.019, closest I could get to 0.05

    Ljung-Box test

data:  Residuals from ARIMA(6,2,9)
Q* = 16.686, df = 7, p-value = 0.01953

Model df: 15.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_UK <- sarima.for(casesUK, n.ahead = 10, 3, 1, 2)

fc_10_UK$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 4248.787 4537.199 4975.123 5093.243 4990.550 4764.910 4639.472 4696.669 4896.270 5103.977

fcd_10_UK <- sarima.for(deathsUK, n.ahead = 10, 6, 2, 9) 

fcd_10_UK$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 237.2505 229.5029 450.6278 506.4696 410.9696 425.0414 204.9640 106.4399 127.9745 236.3315

4. State level time series analysis A. Filter data to 7 key states (MN, CO, MI, NE, PA, SD, TX) B. Replicate ARIMA models for each state C. Forecast 10 days out D. Multivariate ts analysis (VAR and MARSS) for MN

autoplot(fit_casesSD)
autoplot(fit_casesSD)

sarima.for(casesSD_train, n.ahead = 20, 0, 1, 3)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  55.80984  78.14699  86.87687  87.73647  88.59606  89.45565  90.31525  91.17484  92.03444  92.89403
[11]  93.75363  94.61322  95.47281  96.33241  97.19200  98.05160  98.91119  99.77079 100.63038 101.48997

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 14.55914 19.17060 24.74320 32.90091 39.40448 44.97729 49.93196 54.43752 58.59767 62.48145 66.13754
[12] 69.60185 72.90173 76.05856 79.08950 82.00849 84.82709 87.55501 90.20046 92.77051

lines(casesSD_test)

lines(casesSD_test)

checkresiduals(fit_casesSD) # doesn't pass Ljung-Box test

    Ljung-Box test

data:  Residuals from ARIMA(0,1,3)
Q* = 34.763, df = 19, p-value = 0.01491

Model df: 3.   Total lags used: 22

lines(casesSD_test)

checkresiduals(fit_casesSD) # doesn't pass Ljung-Box test

    Ljung-Box test

data:  Residuals from ARIMA(0,1,3)
Q* = 34.763, df = 19, p-value = 0.01491

Model df: 3.   Total lags used: 22

fit_casesSD2 <- arima(casesSD, order = c(2, 1, 5))
checkresiduals(fit_casesSD) # doesn't pass Ljung-Box test

    Ljung-Box test

data:  Residuals from ARIMA(0,1,3)
Q* = 34.763, df = 19, p-value = 0.01491

Model df: 3.   Total lags used: 22

fit_casesSD2 <- arima(casesSD, order = c(2, 1, 5))
checkresiduals(fit_casesSD2) # p-value = 0.018 (closest I could get to 0.05) - not a good model

    Ljung-Box test

data:  Residuals from ARIMA(2,1,5)
Q* = 28.669, df = 15, p-value = 0.01774

Model df: 7.   Total lags used: 22

fit_deathsSD <- auto.arima(deathsSD, stepwise = FALSE, approximation = FALSE)
fit_deathsSD # ARIMA(4, 1, 0), AICc - 217.37

Series: deathsSD 
ARIMA(4,1,0) 

Coefficients:
          ar1      ar2      ar3      ar4
      -0.3889  -0.6590  -0.5299  -0.4510
s.e.   0.0930   0.0832   0.0888   0.1015

sigma^2 estimated as 0.4121:  log likelihood=-103.39
AIC=216.78   AICc=217.37   BIC=230.14

fit_deathsSD # ARIMA(4, 1, 0), AICc - 217.37

Series: deathsSD 
ARIMA(4,1,0) 

Coefficients:
          ar1      ar2      ar3      ar4
      -0.3889  -0.6590  -0.5299  -0.4510
s.e.   0.0930   0.0832   0.0888   0.1015

sigma^2 estimated as 0.4121:  log likelihood=-103.39
AIC=216.78   AICc=217.37   BIC=230.14

autoplot(fit_deathsSD)

autoplot(fit_deathsSD)

sarima.for(deathsSD_train, n.ahead = 20, 4, 1, 0)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 0.6879463 0.8160066 0.7666529 0.2831663 1.2283244 0.9142358 0.7809673 0.9437538 0.6098081 0.9502225
[11] 0.9661577 0.8287993 0.9597614 0.8140346 0.8973339 0.9743423 0.8951964 0.9621354 0.9227950 0.9259231

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 0.3920909 0.4193993 0.4243495 0.4382226 0.4396933 0.4856801 0.5140001 0.5229953 0.5396592 0.5463374
[11] 0.5638882 0.5848724 0.5958470 0.6103421 0.6207092 0.6326704 0.6482213 0.6596789 0.6720797 0.6834272

lines(deathsSD_test)

checkresiduals(fit_deathsSD) # passes Ljung-Box test

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0)
Q* = 27.879, df = 18, p-value = 0.06392

Model df: 4.   Total lags used: 22

lines(deathsSD_test)

checkresiduals(fit_deathsSD) # passes Ljung-Box test

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0)
Q* = 27.879, df = 18, p-value = 0.06392

Model df: 4.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_SD <- sarima.for(casesSD, n.ahead = 10, 2, 1, 5)

fc_10_SD$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 286.1869 314.0617 281.5888 252.0685 231.5716 212.9217 201.0289 189.2848 182.8396 175.4861

fcd_10_SD <- sarima.for(deathsSD, n.ahead = 10, 4, 1, 0) 

fcd_10_SD <- sarima.for(deathsSD, n.ahead = 10, 4, 1, 0) 

fcd_10_SD$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 3.918686 4.102494 2.803120 1.399809 2.378335 3.620602 3.914558 3.172091 2.231339 2.448018

fit_casesSD2 <- arima(casesSD, order = c(2, 1, 5))
checkresiduals(fit_casesSD2) # p-value = 0.018 (closest I could get to 0.05) - not a good model

    Ljung-Box test

data:  Residuals from ARIMA(2,1,5)
Q* = 28.669, df = 15, p-value = 0.01774

Model df: 7.   Total lags used: 22

fit_deathsSD <- auto.arima(deathsSD, stepwise = FALSE, approximation = FALSE)
fit_deathsSD # ARIMA(4, 1, 0), AICc - 217.37

Series: deathsSD 
ARIMA(4,1,0) 

Coefficients:
          ar1      ar2      ar3      ar4
      -0.3889  -0.6590  -0.5299  -0.4510
s.e.   0.0930   0.0832   0.0888   0.1015

sigma^2 estimated as 0.4121:  log likelihood=-103.39
AIC=216.78   AICc=217.37   BIC=230.14

autoplot(fit_deathsSD)

sarima.for(deathsSD_train, n.ahead = 20, 4, 1, 0)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 0.6879463 0.8160066 0.7666529 0.2831663 1.2283244 0.9142358 0.7809673 0.9437538 0.6098081 0.9502225
[11] 0.9661577 0.8287993 0.9597614 0.8140346 0.8973339 0.9743423 0.8951964 0.9621354 0.9227950 0.9259231

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 0.3920909 0.4193993 0.4243495 0.4382226 0.4396933 0.4856801 0.5140001 0.5229953 0.5396592 0.5463374
[11] 0.5638882 0.5848724 0.5958470 0.6103421 0.6207092 0.6326704 0.6482213 0.6596789 0.6720797 0.6834272

lines(deathsSD_test)

checkresiduals(fit_deathsSD) # passes Ljung-Box test

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0)
Q* = 27.879, df = 18, p-value = 0.06392

Model df: 4.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_SD <- sarima.for(casesSD, n.ahead = 10, 2, 1, 5)

fc_10_SD$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 286.1869 314.0617 281.5888 252.0685 231.5716 212.9217 201.0289 189.2848 182.8396 175.4861

fcd_10_SD <- sarima.for(deathsSD, n.ahead = 10, 4, 1, 0) 

fcd_10_SD$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 3.918686 4.102494 2.803120 1.399809 2.378335 3.620602 3.914558 3.172091 2.231339 2.448018

# Repeat for TX
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS, DIFFERENCED DATA SETS
casesTX <- tsTX[, "Cases"]
casesTX_train <- casesTX %>% window(end = c(2020, 120))
casesTX_train <- casesTX %>% window(end = c(2020, 120))
casesTX_test <- casesTX %>% window(end = c(2020, 130))
diff_casesTX <- diff(casesTX)
deathsTX <- tsTX[, "Deaths"]
deathsTX <- tsTX[, "Deaths"]
deathsTX_train <- deathsTX %>% window(end = c(2020, 120))
deathsTX_train <- deathsTX %>% window(end = c(2020, 120))
deathsTX_test <- deathsTX %>% window(end = c(2020, 130))
diff_deathsTX <- diff(deathsTX)
diff_deathsTX <- diff(deathsTX)
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_casesTX <- auto.arima(casesTX, stepwise = FALSE, approximation = FALSE)
fit_casesTX # ARIMA(5, 1, 0), AICc - 1358.83

Series: casesTX 
ARIMA(5,1,0) with drift 

Coefficients:
          ar1      ar2      ar3      ar4      ar5    drift
      -0.6897  -0.0959  -0.3050  -0.6786  -0.3295  10.7237
s.e.   0.0906   0.0910   0.0845   0.0881   0.0907   4.0632

sigma^2 estimated as 17299:  log likelihood=-671.85
AIC=1357.7   AICc=1358.83   BIC=1376.41

# BEST MODELS - use auto.arima models for simplicity and consistency
fit_casesTX <- auto.arima(casesTX, stepwise = FALSE, approximation = FALSE)
fit_casesTX # ARIMA(5, 1, 0), AICc - 1358.83

Series: casesTX 
ARIMA(5,1,0) with drift 

Coefficients:
          ar1      ar2      ar3      ar4      ar5    drift
      -0.6897  -0.0959  -0.3050  -0.6786  -0.3295  10.7237
s.e.   0.0906   0.0910   0.0845   0.0881   0.0907   4.0632

sigma^2 estimated as 17299:  log likelihood=-671.85
AIC=1357.7   AICc=1358.83   BIC=1376.41

autoplot(fit_casesTX)

sarima.for(casesTX_train, n.ahead = 20, 5, 1, 0)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1]  928.6433 1013.9360  614.0468  827.4279  760.7534  858.5194 1044.1534  885.4722  955.4576  842.1049
[11]  822.1472  935.1035  912.2219 1012.1075  980.7167  938.7568  956.7482  918.9534  982.6893 1006.9452

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 116.0575 124.7555 143.6116 146.3068 146.8303 152.0277 153.7902 167.2111 178.8911 181.9872 186.3042
[12] 186.7495 189.3770 194.6034 199.5095 206.9761 210.0200 212.5235 214.8633 217.0844

lines(casesTX_test)

checkresiduals(fit_casesTX) # too low

    Ljung-Box test

data:  Residuals from ARIMA(5,1,0) with drift
Q* = 40.296, df = 16, p-value = 0.0007048

Model df: 6.   Total lags used: 22

fit_casesTX2 <- arima(casesTX, order = c(7, 1, 0))
checkresiduals(fit_casesTX2) # passes the Ljung Box test

    Ljung-Box test

data:  Residuals from ARIMA(7,1,0)
Q* = 22.883, df = 15, p-value = 0.08665

Model df: 7.   Total lags used: 22

fit_deathsTX <- auto.arima(deathsTX, stepwise = FALSE, approximation = FALSE)
fit_deathsTX # ARIMA(5, 1, 0), AICc - 689.06

Series: deathsTX 
ARIMA(5,1,0) with drift 

Coefficients:
          ar1      ar2      ar3      ar4      ar5   drift
      -0.5822  -0.0692  -0.3807  -0.6652  -0.2355  0.3269
s.e.   0.0939   0.0891   0.0792   0.0869   0.0932  0.1877

sigma^2 estimated as 33.07:  log likelihood=-336.96
AIC=687.92   AICc=689.06   BIC=706.63

fit_deathsTX # ARIMA(5, 1, 0), AICc - 689.06

Series: deathsTX 
ARIMA(5,1,0) with drift 

Coefficients:
          ar1      ar2      ar3      ar4      ar5   drift
      -0.5822  -0.0692  -0.3807  -0.6652  -0.2355  0.3269
s.e.   0.0939   0.0891   0.0792   0.0869   0.0932  0.1877

sigma^2 estimated as 33.07:  log likelihood=-336.96
AIC=687.92   AICc=689.06   BIC=706.63

autoplot(fit_deathsTX)

sarima.for(deathsTX_train, n.ahead = 20, 5, 1, 0)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 32.954946 21.172757  6.384597 30.738018 33.830672 41.542596 38.961330 18.641371 22.023326 23.109805
[11] 33.453473 43.970925 34.726711 30.373163 23.161839 24.136606 35.708556 38.441232 39.749723 33.404030

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 4.698068 5.047545 5.337290 5.360786 5.510997 5.847467 6.158233 6.713940 7.103123 7.106117 7.125461
[12] 7.152362 7.362992 7.908017 8.054095 8.130455 8.136884 8.151303 8.363989 8.593711

lines(deathsTX_test)

checkresiduals(fit_deathsTX) # passes

    Ljung-Box test

data:  Residuals from ARIMA(5,1,0) with drift
Q* = 16.177, df = 16, p-value = 0.4407

Model df: 6.   Total lags used: 22

lines(deathsTX_test)

checkresiduals(fit_deathsTX) # passes

    Ljung-Box test

data:  Residuals from ARIMA(5,1,0) with drift
Q* = 16.177, df = 16, p-value = 0.4407

Model df: 6.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_TX <- sarima.for(casesTX, n.ahead = 10, 7, 1, 0)

fc_10_TX$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 1036.677 1062.171 1099.984 1071.121 1248.110 1200.434 1179.260 1168.616 1130.126 1153.625

fcd_10_TX <- sarima.for(deathsTX, n.ahead = 10, 5, 1, 0) 

fcd_10_TX$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 36.43089 34.04652 29.87068 32.38012 35.87844 38.87863 40.23278 38.17785 36.17898 35.10841

V. County level time series analysis A. Filter data to 34 key counties (MN - 9, CO - 8, MI - 4, NE - 4, PA - 3, SD - 2, TX - 4) B. Replicate ARIMA models for each county C. Forecast 10 days out

# BEST MODELS - use auto.arima models for simplicity and consistency
fit_Tarrant_cases <- auto.arima(Tarrant_cases, stepwise = FALSE, approximation = FALSE)
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_Tarrant_cases <- auto.arima(Tarrant_cases, stepwise = FALSE, approximation = FALSE)
fit_Tarrant_cases # ARIMA(3, 1, 2), AICc - 1076.92

Series: Tarrant_cases 
ARIMA(3,1,2) with drift 

Coefficients:
         ar1     ar2      ar3      ma1     ma2   drift
      0.6389  0.1859  -0.2497  -1.8684  0.9229  0.9503
s.e.  0.1086  0.1168   0.1208   0.0662  0.0556  0.4227

sigma^2 estimated as 1219:  log likelihood=-530.89
AIC=1075.78   AICc=1076.92   BIC=1094.49

autoplot(fit_Tarrant_cases)

autoplot(fit_Tarrant_cases)

sarima.for(Tarrant_cases_train, n.ahead = 20, 3, 1, 2)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 89.82678 71.49481 90.46510 77.39054 89.05233 82.29664 88.62381 85.85879 89.23810 88.53458 90.48947
[12] 90.73914 92.06385 92.72083 93.78557 94.60540 95.56952 96.45040 97.37808 98.28026

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 27.20873 28.37068 28.69945 29.90789 29.92150 30.68212 30.71986 31.17960 31.27613 31.58981 31.73717
[12] 31.98524 32.16164 32.38046 32.56967 32.77506 32.96846 33.16715 33.36083 33.55561

autoplot(fit_Tarrant_cases)

sarima.for(Tarrant_cases_train, n.ahead = 20, 3, 1, 2)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 89.82678 71.49481 90.46510 77.39054 89.05233 82.29664 88.62381 85.85879 89.23810 88.53458 90.48947
[12] 90.73914 92.06385 92.72083 93.78557 94.60540 95.56952 96.45040 97.37808 98.28026

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 27.20873 28.37068 28.69945 29.90789 29.92150 30.68212 30.71986 31.17960 31.27613 31.58981 31.73717
[12] 31.98524 32.16164 32.38046 32.56967 32.77506 32.96846 33.16715 33.36083 33.55561

lines(Tarrant_cases_test)

lines(Tarrant_cases_test)

checkresiduals(fit_Tarrant_cases) # does not pass Ljung Box test

    Ljung-Box test

data:  Residuals from ARIMA(3,1,2) with drift
Q* = 31.754, df = 16, p-value = 0.01076

Model df: 6.   Total lags used: 22

fit_Tarrant_cases2 <- arima(Tarrant_cases, order = c(3, 1, 8))
checkresiduals(fit_Tarrant_cases) # does not pass Ljung Box test

    Ljung-Box test

data:  Residuals from ARIMA(3,1,2) with drift
Q* = 31.754, df = 16, p-value = 0.01076

Model df: 6.   Total lags used: 22

fit_Tarrant_cases2 <- arima(Tarrant_cases, order = c(3, 1, 8))
checkresiduals(fit_Tarrant_cases2) # passes

    Ljung-Box test

data:  Residuals from ARIMA(3,1,8)
Q* = 18.296, df = 11, p-value = 0.07496

Model df: 11.   Total lags used: 22

fit_Tarrant_deaths <- auto.arima(Tarrant_deaths, stepwise = FALSE, approximation = FALSE)
fit_Tarrant_deaths # ARIMA(2, 1, 3), AICc - 375.33

Series: Tarrant_deaths 
ARIMA(2,1,3) 

Coefficients:
          ar1     ar2      ma1      ma2     ma3
      -0.2626  0.2820  -0.6236  -0.8434  0.7056
s.e.   0.1527  0.1439   0.1137   0.0602  0.1032

sigma^2 estimated as 1.764:  log likelihood=-181.24
AIC=374.49   AICc=375.33   BIC=390.52

fit_Tarrant_deaths <- auto.arima(Tarrant_deaths, stepwise = FALSE, approximation = FALSE)
fit_Tarrant_deaths # ARIMA(2, 1, 3), AICc - 375.33

Series: Tarrant_deaths 
ARIMA(2,1,3) 

Coefficients:
          ar1     ar2      ma1      ma2     ma3
      -0.2626  0.2820  -0.6236  -0.8434  0.7056
s.e.   0.1527  0.1439   0.1137   0.0602  0.1032

sigma^2 estimated as 1.764:  log likelihood=-181.24
AIC=374.49   AICc=375.33   BIC=390.52

autoplot(fit_Tarrant_deaths)

fit_Tarrant_deaths # ARIMA(2, 1, 3), AICc - 375.33

Series: Tarrant_deaths 
ARIMA(2,1,3) 

Coefficients:
          ar1     ar2      ma1      ma2     ma3
      -0.2626  0.2820  -0.6236  -0.8434  0.7056
s.e.   0.1527  0.1439   0.1137   0.0602  0.1032

sigma^2 estimated as 1.764:  log likelihood=-181.24
AIC=374.49   AICc=375.33   BIC=390.52

autoplot(fit_Tarrant_deaths)

sarima.for(Tarrant_deaths_train, n.ahead = 20, 2, 1, 3)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 1.437549 2.375279 1.894045 2.274475 2.086029 2.250629 2.185897 2.265192 2.251753 2.296933 2.304551
[12] 2.335998 2.352216 2.378113 2.397833 2.421484 2.442628 2.465368 2.487091 2.509461

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 0.9729671 0.9738641 0.9928969 0.9974776 1.0140056 1.0234201 1.0377454 1.0489840 1.0620761 1.0738322
[11] 1.0862361 1.0980271 1.1099984 1.1216529 1.1333085 1.1447699 1.1561657 1.1674202 1.1785863 1.1896355

lines(Tarrant_deaths_test)

sarima.for(Tarrant_deaths_train, n.ahead = 20, 2, 1, 3)

$pred
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 1.437549 2.375279 1.894045 2.274475 2.086029 2.250629 2.185897 2.265192 2.251753 2.296933 2.304551
[12] 2.335998 2.352216 2.378113 2.397833 2.421484 2.442628 2.465368 2.487091 2.509461

$se
Time Series:
Start = c(2020, 121) 
End = c(2020, 140) 
Frequency = 365 
 [1] 0.9729671 0.9738641 0.9928969 0.9974776 1.0140056 1.0234201 1.0377454 1.0489840 1.0620761 1.0738322
[11] 1.0862361 1.0980271 1.1099984 1.1216529 1.1333085 1.1447699 1.1561657 1.1674202 1.1785863 1.1896355

lines(Tarrant_deaths_test)

checkresiduals(fit_Tarrant_deaths) # does not pass

    Ljung-Box test

data:  Residuals from ARIMA(2,1,3)
Q* = 42.544, df = 17, p-value = 0.0005594

Model df: 5.   Total lags used: 22

fit_Tarrant_deaths2 <- arima(Tarrant_deaths, order = c(2, 2, 6))
checkresiduals(fit_Tarrant_deaths2) # passes

    Ljung-Box test

data:  Residuals from ARIMA(2,2,6)
Q* = 22.922, df = 14, p-value = 0.06155

Model df: 8.   Total lags used: 22

fit_Tarrant_deaths2 <- arima(Tarrant_deaths, order = c(2, 2, 6))
checkresiduals(fit_Tarrant_deaths2) # passes

    Ljung-Box test

data:  Residuals from ARIMA(2,2,6)
Q* = 22.922, df = 14, p-value = 0.06155

Model df: 8.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_Tarrant <- sarima.for(Tarrant_cases, n.ahead = 10, 3, 1, 8)

# Use best models to forecast further ahead
fc_10_Tarrant <- sarima.for(Tarrant_cases, n.ahead = 10, 3, 1, 8)
fc_10_Tarrant$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 116.19146  89.58502 150.79515  32.90965 177.56836  81.34822 119.76881 130.77723  96.63236 138.87275

lines(Tarrant_deaths_test)

checkresiduals(fit_Tarrant_deaths) # does not pass

    Ljung-Box test

data:  Residuals from ARIMA(2,1,3)
Q* = 42.544, df = 17, p-value = 0.0005594

Model df: 5.   Total lags used: 22

fit_Tarrant_deaths2 <- arima(Tarrant_deaths, order = c(2, 2, 6))
checkresiduals(fit_Tarrant_deaths2) # passes

    Ljung-Box test

data:  Residuals from ARIMA(2,2,6)
Q* = 22.922, df = 14, p-value = 0.06155

Model df: 8.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_Tarrant <- sarima.for(Tarrant_cases, n.ahead = 10, 3, 1, 8)

fc_10_Tarrant$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 116.19146  89.58502 150.79515  32.90965 177.56836  81.34822 119.76881 130.77723  96.63236 138.87275

fcd_10_Tarrant <- sarima.for(Tarrant_deaths, n.ahead = 10, 2, 2, 6) 
lines(Tarrant_deaths_test)

checkresiduals(fit_Tarrant_deaths) # does not pass

    Ljung-Box test

data:  Residuals from ARIMA(2,1,3)
Q* = 42.544, df = 17, p-value = 0.0005594

Model df: 5.   Total lags used: 22

fit_Tarrant_deaths2 <- arima(Tarrant_deaths, order = c(2, 2, 6))
checkresiduals(fit_Tarrant_deaths2) # passes

    Ljung-Box test

data:  Residuals from ARIMA(2,2,6)
Q* = 22.922, df = 14, p-value = 0.06155

Model df: 8.   Total lags used: 22

# Use best models to forecast further ahead
fc_10_Tarrant <- sarima.for(Tarrant_cases, n.ahead = 10, 3, 1, 8)

fc_10_Tarrant$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 116.19146  89.58502 150.79515  32.90965 177.56836  81.34822 119.76881 130.77723  96.63236 138.87275

fcd_10_Tarrant <- sarima.for(Tarrant_deaths, n.ahead = 10, 2, 2, 6) 

fcd_10_Tarrant$pred

Time Series:
Start = c(2020, 131) 
End = c(2020, 140) 
Frequency = 365 
 [1] 3.126920 5.342878 3.640292 3.691021 5.676126 2.693883 5.921924 3.462516 4.858318 4.844965

V. Linear Regression on MN counties to account for constants A. Load current data and add constant variables B. Reshape data into long and wide forms C. Build linear regression models for cases and deaths

# load in MN county level data from usafacts (wide form)
# http://www.cookbook-r.com/Manipulating_data/Converting_data_between_wide_and_long_format/
MN_0613 <- MN_0613 %>% filter(State == "MN") %>% filter(FIPS != "0") # remove unallocated cases

Error: Problem with `filter()` input `..1`.
x object 'State' not found
ℹ Input `..1` is `State == "MN"`.

6. Vector Autoregression (VAR) A. Perform VAR on MN counties (9) and select counties (3) B. Perform VAR on US states (7) and select states(3) C. Forecast 12 days out on all models

7. Final model check A. Load current data through June 24th (COVID cases and deaths) B. Check predictions from 3 key county models against current data (ARIMA, VAR, Regression) C. Check predictions for state of MN against current data (ARIMA, VAR)

# ADD MORE RECENT DATA TO US, MN, SELECT COUNTIES (FROM USAFACTS)
# RERUN TS FORECASTS USING UPDATED DATA AND MEASURE ACCURACY, PREDICT NEXT 10 DAYS

# data through 06/12/2020
# https://usafacts.org/visualizations/coronavirus-covid-19-spread-map/

# steps to combine data
# filter cases to single county
MN0612c <- cases0612 %>% filter(State == "MN")
Hennepin0612c <- MN0612c %>% filter(County.Name == "Hennepin County")
Hennepin0612c <- Hennepin0612c[, 5:148]
dim(Hennepin0612c)
class(Hennepin0612c)
# invert columns and rows
x <- as.data.frame(t(Hennepin0612c))
dim(x)
class(x)
colnames(x) <- "V1"
str(x)
start <- x[1, 1]
start
tail(x)
# add col that shows difference between values (USAfacts data is cumulative)
x$diff <- c(start, diff(x$V1))
x
colnames(x) <- c("Cum", "Cases")
# repeat for deaths
MN0612d <- deaths0612 %>% filter(State == "MN")
Hennepin0612d <- MN0612d %>% filter(County.Name == "Hennepin County")
Hennepin0612d <- Hennepin0612d[, 5:148]
str(Hennepin0612d)
dim(Hennepin0612d)
# invert columns and rows
y <- as.data.frame(t(Hennepin0612d))
dim(y)
colnames(y) <- "V1"
str(y)
start <- y[1, 1]
start
tail(y)
# add col that shows difference between values (USAfacts data is cumulative)
y$diff <- c(start, diff(y$V1))
head(y)
tail(y)
colnames(y) <- c("Cum", "Deaths")
# combine Cases and Deaths
HennMN0612 <- cbind(x$Cases, y$Deaths)
colnames(HennMN0612) <- c("Cases", "Deaths")
HennMN0612
str(HennMN0612)
tail(HennMN0612)
HennMNts <- ts(HennMN0612, start = c(2020, 22), frequency = 365)
str(HennMNts)
head(HennMNts)
HennMNts
autoplot(HennMNts, main = "COVID-19 Cases and Deaths - Minnesota")


# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS
autoplot(HennMNts)
Henn_cases <- HennMNts[, "Cases"]
Henn_cases_train <- Henn_cases %>% window(end = c(2020, 124))
Henn_cases_test <- Henn_cases %>% window(end = c(2020, 144))
Henn_deaths <- HennMNts[, "Deaths"]
Henn_deaths_train <- Henn_deaths %>% window(end = c(2020, 124))
Henn_deaths_test <- Henn_deaths %>% window(end = c(2020, 144))
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_Henn_cases <- auto.arima(Henn_cases, stepwise = FALSE, approximation = FALSE)
fit_Henn_cases # ARIMA(0, 1, 1), AICc - 1603.27
autoplot(fit_Henn_cases)
sarima.for(Henn_cases_train, n.ahead = 20, 0, 1, 1)
lines(Henn_cases_test)
checkresiduals(fit_Henn_cases) # passes
fit_Henn_deaths <- auto.arima(Henn_deaths, stepwise = FALSE, approximation = FALSE)
fit_Henn_deaths # ARIMA(3, 1, 2), AICc - 725.33
autoplot(fit_Henn_deaths)
sarima.for(Henn_deaths_train, n.ahead = 20, 3, 1, 2)
lines(Henn_deaths_test)
checkresiduals(fit_Henn_deaths) # p-value too low, does not pass Ljung-Box test
fit_Henn_deaths2 <- arima(Henn_deaths, order = c(6, 1, 1))
sarima.for(Henn_deaths_train, n.ahead = 20, 6, 1, 1)
lines(Henn_deaths_test)
checkresiduals(fit_Henn_deaths2) # passes
# Use best models to forecast further ahead
fc_10_Henn <- sarima.for(Henn_cases, n.ahead = 10, 0, 1, 1)
fc_10_Henn$pred
fcd_10_Henn <- sarima.for(Henn_deaths, n.ahead = 10, 6, 1, 1) 
fcd_10_Henn$pred

# REPEAT PROCESS
# steps to combine data
# filter cases to single county
MN0612c <- cases0612 %>% filter(State == "MN")
Stearns0612c <- MN0612c %>% filter(County.Name == "Stearns County")
Stearns0612c <- Stearns0612c[, 5:148]
dim(Stearns0612c)
# invert columns and rows
x <- as.data.frame(t(Stearns0612c))
colnames(x) <- "V1"
start <- x[1, 1]
start
# add col that shows difference between values (USAfacts data is cumulative)
x$diff <- c(start, diff(x$V1))
colnames(x) <- c("Cum", "Cases")
# repeat for deaths
MN0612d <- deaths0612 %>% filter(State == "MN")
Stearns0612d <- MN0612d %>% filter(County.Name == "Stearns County")
Stearns0612d <- Stearns0612d[, 5:148]
dim(Stearns0612d)
# invert columns and rows
y <- as.data.frame(t(Stearns0612d))
colnames(y) <- "V1"
start <- y[1, 1]
start
# add col that shows difference between values (USAfacts data is cumulative)
y$diff <- c(start, diff(y$V1))
colnames(y) <- c("Cum", "Deaths")
# combine Cases and Deaths
StearnsMN0612 <- cbind(x$Cases, y$Deaths)
colnames(StearnsMN0612) <- c("Cases", "Deaths")
StearnsMNts <- ts(StearnsMN0612, start = c(2020, 22), frequency = 365)
autoplot(StearnsMNts, main = "COVID-19 Cases and Deaths - MN, Stearns County")
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS
autoplot(StearnsMNts)
Stearns_cases <- StearnsMNts[, "Cases"]
Stearns_cases_train <- Stearns_cases %>% window(end = c(2020, 124))
Stearns_cases_test <- Stearns_cases %>% window(end = c(2020, 144))
Stearns_deaths <- StearnsMNts[, "Deaths"]
Stearns_deaths_train <- Stearns_deaths %>% window(end = c(2020, 124))
Stearns_deaths_test <- Stearns_deaths %>% window(end = c(2020, 144))
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_Stearns_cases <- auto.arima(Stearns_cases, stepwise = FALSE, approximation = FALSE)
fit_Stearns_cases # ARIMA(4, 1, 0), AICc - 1197.98
autoplot(fit_Stearns_cases)
sarima.for(Stearns_cases_train, n.ahead = 20, 4, 1, 0)
lines(Stearns_cases_test)
checkresiduals(fit_Stearns_cases) # p-value too low
fit_Stearns_cases2 <- arima(Stearns_cases, order = c(4, 1, 2))
checkresiduals(fit_Stearns_cases2) # passes
fit_Stearns_deaths <- auto.arima(Stearns_deaths, stepwise = FALSE, approximation = FALSE)
fit_Stearns_deaths # ARIMA(0, 1, 5), AICc - 94.28
autoplot(fit_Stearns_deaths)
sarima.for(Stearns_deaths_train, n.ahead = 20, 0, 1, 5)
checkresiduals(fit_Stearns_deaths) # p-value too low, does not pass Ljung-Box test
fit_Stearns_deaths2 <- arima(Stearns_deaths, order = c(0, 1, 6))
checkresiduals(fit_Stearns_deaths2) # passes
# Use best models to forecast further ahead
fc_10_Stearns <- sarima.for(Stearns_cases, n.ahead = 10, 4, 1, 2)
fc_10_Stearns$pred
fcd_10_Stearns <- sarima.for(Stearns_deaths, n.ahead = 10, 0, 1, 6) 
fcd_10_Stearns$pred

# REPEAT PROCESS
# steps to combine data
# filter cases to single county
MN0612c <- cases0612 %>% filter(State == "MN")
Kandiyohi0612c <- MN0612c %>% filter(County.Name == "Kandiyohi County")
Kandiyohi0612c <- Kandiyohi0612c[, 5:148]
dim(Kandiyohi0612c)
# invert columns and rows
x <- as.data.frame(t(Kandiyohi0612c))
colnames(x) <- "V1"
start <- x[1, 1]
start
# add col that shows difference between values (USAfacts data is cumulative)
x$diff <- c(start, diff(x$V1))
colnames(x) <- c("Cum", "Cases")
# repeat for deaths
MN0612d <- deaths0612 %>% filter(State == "MN")
Kandiyohi0612d <- MN0612d %>% filter(County.Name == "Kandiyohi County")
Kandiyohi0612d <- Kandiyohi0612d[, 5:148]
dim(Kandiyohi0612d)
# invert columns and rows
y <- as.data.frame(t(Kandiyohi0612d))
colnames(y) <- "V1"
start <- y[1, 1]
start
# add col that shows difference between values (USAfacts data is cumulative)
y$diff <- c(start, diff(y$V1))
colnames(y) <- c("Cum", "Deaths")
# combine Cases and Deaths
KandiyohiMN0612 <- cbind(x$Cases, y$Deaths)
colnames(KandiyohiMN0612) <- c("Cases", "Deaths")
KandiyohiMNts <- ts(KandiyohiMN0612, start = c(2020, 22), frequency = 365)
autoplot(KandiyohiMNts, main = "COVID-19 Cases and Deaths - MN, Kandiyohi County")
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS
autoplot(KandiyohiMNts)
Kandiyohi_cases <- KandiyohiMNts[, "Cases"]
Kandiyohi_cases_train <- Kandiyohi_cases %>% window(end = c(2020, 124))
Kandiyohi_cases_test <- Kandiyohi_cases %>% window(end = c(2020, 144))
Kandiyohi_deaths <- KandiyohiMNts[, "Deaths"]
Kandiyohi_deaths_train <- Kandiyohi_deaths %>% window(end = c(2020, 124))
Kandiyohi_deaths_test <- Kandiyohi_deaths %>% window(end = c(2020, 144))
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_Kandiyohi_cases <- auto.arima(Kandiyohi_cases, stepwise = FALSE, approximation = FALSE)
fit_Kandiyohi_cases # ARIMA(2, 1, 3), AICc - 914.05
autoplot(fit_Kandiyohi_cases)
sarima.for(Kandiyohi_cases_train, n.ahead = 20, 2, 1, 3)
lines(Kandiyohi_cases_test)
checkresiduals(fit_Kandiyohi_cases) # passes
fit_Kandiyohi_deaths <- auto.arima(Kandiyohi_deaths, stepwise = FALSE, approximation = FALSE)
fit_Kandiyohi_deaths # ARIMA(0, 0, 0), AICc - -304.97 White noise model
autoplot(fit_Kandiyohi_deaths)
sarima.for(Kandiyohi_deaths_train, n.ahead = 20, 0, 0, 0)
checkresiduals(fit_Kandiyohi_deaths) # passes Ljung-Box test
# Use best models to forecast further ahead
fc_10_Kandiyohi <- sarima.for(Kandiyohi_cases, n.ahead = 10, 2, 1, 3)
fc_10_Kandiyohi$pred
fcd_10_Kandiyohi <- sarima.for(Kandiyohi_deaths, n.ahead = 10, 0, 0, 0) 
fcd_10_Kandiyohi$pred

# REPEAT for MN Data (full state)
MN0612c <- cases0612 %>% filter(State == "MN")
MN0612c <- MN0612c[, 5:148]
dim(MN0612c)
total <- colSums(MN0612c)
length(total)
total <- as.integer(total)
MN0612c <- rbind(MN0612c, total)
class(MN0612c)
MN0612c <- MN0612c[89, ]
MN0612c
# invert columns and rows
x <- as.data.frame(t(MN0612c))
class(x)
str(x)
colnames(x) <- "V1"
str(x)
start <- x[1, 1]
start
# add col that shows difference between values (USAfacts data is cumulative)
x$diff <- c(start, diff(x$V1))
colnames(x) <- c("Cum", "Cases")
tail(x)
# repeat for deaths (full state)
MN0612d <- deaths0612 %>% filter(State == "MN")
MN0612d <- MN0612d[, 5:148]
dim(MN0612d)
total <- colSums(MN0612d)
length(total)
total <- as.integer(total)
MN0612d <- rbind(MN0612d, total)
class(MN0612d)
MN0612d <- MN0612d[89, ]
MN0612d
# invert columns and rows
y <- as.data.frame(t(MN0612d))
colnames(y) <- "V1"
str(y)
start <- y[1, 1]
start
# add col that shows difference between values (USAfacts data is cumulative)
y$diff <- c(start, diff(y$V1))
colnames(y) <- c("Cum", "Deaths")
tail(y)
# combine Cases and Deaths
MN0612 <- cbind(x$Cases, y$Deaths)
colnames(MN0612) <- c("Cases", "Deaths")
MNts <- ts(MN0612, start = c(2020, 22), frequency = 365)
autoplot(MNts, main = "COVID-19 Cases and Deaths - MN")
# CREATE TEST/TRAIN SETS FOR CASES AND DEATHS
autoplot(MNts)
MN_cases <- MNts[, "Cases"]
MN_cases_train <- MN_cases %>% window(end = c(2020, 124))
MN_cases_test <- MN_cases %>% window(end = c(2020, 144))
MN_deaths <- MNts[, "Deaths"]
MN_deaths_train <- MN_deaths %>% window(end = c(2020, 124))
MN_deaths_test <- MN_deaths %>% window(end = c(2020, 144))
# BEST MODELS - use auto.arima models for simplicity and consistency
fit_MN_cases <- auto.arima(MN_cases, stepwise = FALSE, approximation = FALSE)
fit_MN_cases # ARIMA(3, 1, 2), AICc - 1716.48
autoplot(fit_MN_cases)
sarima.for(MN_cases_train, n.ahead = 20, 3, 1, 2)
lines(MN_cases_test)
checkresiduals(fit_MN_cases) # passes
fit_MN_deaths <- auto.arima(MN_deaths, stepwise = FALSE, approximation = FALSE)
fit_MN_deaths # ARIMA(3, 1, 2), AICc - 836.77
autoplot(fit_MN_deaths)
sarima.for(MN_deaths_train, n.ahead = 20, 3, 1, 2)
lines(MN_deaths_test)
checkresiduals(fit_MN_deaths) # does not pass
fit_MN_deaths2 <- arima(MN_deaths, order = c(4, 1, 3))
checkresiduals(fit_MN_deaths2) # passes
sarima.for(MN_deaths_train, n.ahead = 20, 4, 1, 3)
lines(MN_deaths_test)
# Use best models to forecast further ahead
fc_10_MN <- sarima.for(MN_cases, n.ahead = 10, 3, 1, 2)
fc_10_MN$pred
fcd_10_MN <- sarima.for(MN_deaths, n.ahead = 10, 4, 1, 3) 
fcd_10_MN$pred